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ABSTRACT 

Many analyses and parameter estimations undertaken in astronomy require a large set (> 10 5 ) of 
non-analytical, theoretical spectra, each of these defined by multiple parameters. We describe the 
construction of an Y-dimensional grid which is suitable for generating such spectra. The theoretical 
spectra are designed to correspond to a targeted parameter grid but otherwise to random positions 
in the parameter space, and they are interpolated on-the-fly through a pre-calculated grid of spectra. 
The initial grid is designed to be relatively low in parameter resolution and small in occupied hard 
disk space and therefore can be updated efficiently when a new model is desired. In a pilot study of 
stellar population synthesis of galaxies, the mean square errors on the estimated parameters are found 
to decrease with the targeted grid resolution. This scheme of generating a large model grid is general 
for other areas of studies, particularly if they are based on multi-dimensional parameter space and are 
focused on contrasting model differences. 

Subject headings: techniques: spectroscopic — methods: data analysis 
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1. INTRODUCTION 

Various analyses and parameter estimations performed 
in growing number of applications in astronomy re- 
quire a large set of non-analytical, theoretical spec- 
tra, where each spectrum in-principle can be defined by 
multiple free parameters. In stellar p opulation synthe- 
sis us ing the Bayesian approach (e.g., [K auffman n et all 
2003; GallazzIiranE005; Salim et al. 20017 as many as 
10 4 — I0 5 theoretical spectra and the derived line indices 
are used in order to cover a wide range of star formation 
histories from early- to late-type galaxies. With about 5 
or more free parameters, each theoretical spectrum is a 
single stellar burst defined by metallicity and age, super- 
imposed at a given mass fraction on a spectral compo- 
nent with exponentially decreasing star formation rate 
characterized by the age of the oldest stars, e-folding 
time of star formation, and both spectral components 
could be dust-attenuated. Similarly, 1 5 — I0 6 theoret- 
ical spectra are con sidered in parameter e stimation on 
stellar spectra (e.g.. iRobitaille et al.l I2007L who consid- 
ered a model defined by 14 parameters), and on H II 
regions (e.g.. iMorissetl [20091 who studied a model de- 
fined by 15 parameters). One approach to prepare and 
manage a non-analytical model is to store all of the pre- 
calculated spectra, fixed at both the parameter choice 
and the parameter resolution. As models are becoming 
more sophisticated and are growing in varieties, however, 
updating and storing 10 5 or more spectra may not be the 
most flexible approach for many astronomers. To han- 
dle large sets of theoretical spectra at multi-dimensional 
parameter space is therefore a question that cannot be 
neglected. 

We consider here current integrated stellar population 
models as a case study. These models, together with 
parameter estimation techniques or other analyses, 
were shown by many authors to be invaluable for 
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deriving composition and star formation rate/history 
of different types of galaxies through their observed 
spectra, regardless whether or not a large model 
grid is used (e.g.. iTrager et al.l l2000t iReichardt et al.l 
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20071: l Asari et al.l 120071: iToieiro et al.l 120071 : iChen et al 
2009; Richard s et al.l 120091 ) . Firstly, there is a variety 
of spectral synthesis computational programs which 
use different flavors in initial stellar mass function, 
stellar types, and/or stellar ev o lutionary tracks (e.g. , 



Fioc fe Rocca-Volrnerangel 119971: ILeitherer et al.l 1199 
Bruzual fc CharlotJ |2003t iMarastonl 120051 : iCoelho et al 
20071) . The variations in input ingredients exist not only 
among the models, but also within a single model in the 
form of additional freedom in parameter choice. Sec- 
ondly, increasingly more parameters are being included 
in the analyses, meaning that the parameter space defin- 
ing the theoretical sp ectra is getting larger. So much so 
that iLee et all (|2009f) have recently constructed a stellar 
population model to the level of individual element 
abundances, with the goal to understand the effect of 
each and every element to the integrated spectrum of a 
stellar population. However, there are very few studies 
in the astronomy literature which address the case where 
a multi-dimensional parameter space defines a model. 
Two of the fundamental questions are: (1) how do we 
obtain, or interpolate, a theoretical spectrum in-between 
the parameter values, that is originally unavailable? (2) 
what is the best way for such an interpolation, taking 
into account of the noise in the data, and the possible 
non-linear dependence (jVanderplas fc Connolly! 120091 ) 
between the theoretical spectrum and its underlying 
physical parameters? One of the goals of this work is to 
address the first question. 
Here we describe a novel approach which increases the 
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flexibility of handling a large set (> 10 5 ) of theoreti- 
cal spectra, and is general for TV-dimension (NT)) pa- 
rameter space. The approach adopts a hypercube, an 
NT) analog of a cube of length unity for each side (e.g., 
IAnthony| [l987'). to represent the parameter space under- 
lying each and every theoretical spectrum. The spec- 
trum at any intermediate parameter point is generated 
on-the-ny 2 through multi-linear interpolation, upon an 
initial model grid that is designed to be lower in param- 
eter resolution and therefore can be updated efficiently. 
We examine a pilot case in stellar population synthesis 
of galaxies, to show that the approach is applicable to 
studies in astronomy, and well manageable by a typi- 
cal personal computer - a computer with Intel(R) Pen- 
tium(R)D 3.20 GHz CPU and 3.19 GHz RAM is used in 
this work. 

All of the spectra considered in this work are expressed 
in vacuum wavelength, and are re-sampled to within the 
optical wavelength range 3450 — 8350 A at a resolution 
of 1 A per wavelength bin. The Cartesian coordinate 
system is used throughout this work in representing the 
parameter space. 

2. METHOD: AN TV-DIMENSIONAL GRID FOR MODEL 
PARAMETERS 

Several desirable characteristics are identified in the 
setup for generating a large model grid: (1) generality 
for ND parameter space; (2) the theoretical spectra de- 
fined by the targeted parameters are generated on-the- 
fly, which is essential for improving flexibility and saving 
disk space, especially in arbitrarily complex models; (3) 
targeted parameters can be specified at arbitrary posi- 
tion in the parameter space at run time, with the goal 
to obtain in the analysis/parameter estimation an arbi- 
trary computational resolution which is limited only by 
the parameter resolution in the initial model grid. 

To fulfill the above criteria, we generate a small- to 
moderate-size sample of theoretical spectra to start with. 
The underlying parameters of this initial sample is in the 
form of a Cartesian grid. For each parameter the size of 
the bins is allowed to be uneven, in which case the grid 
will be rectangular instead of square in the 2D analogy. 
For convenience we called this sample the initial grid of 
theoretical spectra, the purpose of which is to form a 
leverage for generating a larger sample of spectra, the 
targeted sample, for which the defining parameters can 
be located anywhere in-between the initial grid points. 
If this sample is also decided to be described by a pa- 
rameter grid instead of random points within the initial 
grid, we called that the targeted grid of spectra. The 
advantage of having a targeted grid of sample is their of- 
fering an intuitive setting for performing necessary inte- 
grations over a parameter space in the applied parameter 
estimation te chnique, an appro a ch also used by various 
authors (e.g., IConti et alJ 120031: iKauffmann et all 120031 : 
IGallazzi et atll2005l : IWelikala et al.H2008fi ~ 

The theoretical spectrum at any intermediate grid 
point is calculated by first locating the corresponding 
neighboring spectra from the initial grid. The neigh- 
boring parameter grid points are then mapped into the 

2 We refer the on-the-fly ability to be the generation of theoreti- 
cal spectra at run time as required by the actual computation, e.g., 
in parameter estimation. 
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Fig. 1. — The time for constructing a hypercube of dimension 
TV, which is used to represent the number of parameters in the 
model. For dust-attenuated integrated stellar spectra defined by 
an exponentially decreasing star formation rate, TV = 4 {age, Z, r e , 
E(B — V)), the time is less than 1 milli-second. For an additional 
dust-attenuated stellar burst, TV = 8 (age, Z, E(B — V) and the 
mass fraction of the burst), the time is less than 1 second. Each 
error bar shows the ± 1-sigma sample scatter of the required time, 
after 20 realizations. The dotted lines are for eye-guiding only. 

corners of a hypercube of dimension TV, where N is the 
number of parameters defining both the initial and tar- 
geted grids. The details of constructing a hypercube is 
described in i j2.ll Multi-linear interpolation (described 
in H2.2p is next used to derive the concerned theoreti- 
cal value at the targeted parameter point. The value is, 
in the current context, the flux density in a given wave- 
length bin. Therefore, the above procedure is repeated 
for each and every wavelength bin of interest. The whole 
procedure of model generation is joined seamlessly with 
the actual computational routines for parameter estima- 
tion or other analyses, with the only input being the 
initial model grid. 

2.1. Hypercube 

A hypercube is an NT) cube with each side equals unity. 
In 2 and 3 dimensions, a hypercube is hence a square 
and a cube, respectively. This geometry makes a hy- 
percube to be naturally suitable for describing an TVD 
Cartesian space, such as a grid of model parameters. For 
our purpose, the dimension of the parameter space is 
defined to be the number of parameters that fully spec- 
ifies a spectrum. For example, in a simple stellar pop- 
ulation defined by stellar age and stellar metallicity, TV 
= 2. In practice, a hypercube is generated by specify- 
ing all the corners using vectors, or arrays in the actual 
computational routines. For example, in 2D the arrays 
are {0, 0}, {0, 1}, {1, 0} and {1, 1}, where any coordinate 
value is defined to be either or 1. It is easy to show by 
deduction that there are in total 2 N corners in an TVD 
hypercube. 

We use a random rendering of points, together with 
a book-keeping approach, to create each and every cor- 
ner of the hypercube exactly once. After some algebraic 
derivations in the low-dimensional cases, we found that 
the number of combinations - or the number of corners 
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Fig. 2. — An illustration of the weighting scheme in the multi- 
linear interpolation, in a 2D hypercube (i.e., a square) which en- 
closes the targeted point (x,y), under the Cartesian coordinate 
system. The weight of point (xo,yo) to the targeted point value is 
the area of the stippled square, \&xq ■ 5yo\. When the point (x, y) is 
moved toward (£2,1/2), the stippled area shrinks, to the point that 
when (x, y) coincides with (x2, y2), &xq ■ Syo = 0, or the weight on 
the targeted point value due to point (xq, yo) is zero. Applying to 
theoretical simple stellar population, for example, the x direction 
can be age, the y direction be Z. 

in the hypercube of dimension N, for k occurrences 
of l's in the array specifying the corner is given by the 
binomial coefficient 



M(N-k)\ 



(1) 



In 3D, there is 1 (= 3C0) combination of no 1: {0,0,0}, 
3 (= 3C1) combinations of one 1: {0,0, 1}, {0, 1,0} and 
{1,0,0}, 3 (= 3C2) combinations of two l's: {0,1,1}, 
{1, 1, 0} and {1, 0, 1}, and 1 (— 3C3) combination of three 
l's: {1,1,1}. Other approaches certainly can be used to 
accomplish the same purpose, as long as a hypercube is 
automatically created given only its dimensionality, N. 

The time required to construct a hypercube of a given 
dimension is shown in Figure [T] For dust-attenuated in- 
tegrated stellar spectra defined by an exponentially de- 
creasing star formation rate, N = 4 {age, Z ', r e , E(B—V), 
at fixed initial stellar mass function and other parame- 
ters), the time is less than 1 milli-second. For an extra 
dust-attenuated stellar burst which is added at a cer- 
tain flux fraction, N — 8 (age, Z, E(B — V) and the 
burst mass relative to the total galaxy mass), the time 
is less than 1 second. The construction of the hypercube 
therefore should not constitute as bottleneck in typical 
analyses or parameter estimation problems, which take 
minutes to hours of computation. 

2.2. Multi-linear interpolation at intermediate grid point 

To interpolate a spectrum at an intermediate parame- 
ter point in the initial grid, called the targeted point, we 
first perform a search to locate all of points, or the neigh- 



boring points, from the initial grid that encloses the tar- 
geted point. This search is performed on the parameter- 
to-parameter basis, for each parameter a pair of neigh- 
boring points is obtained. As a result, the searching time 
scales as N-N gr id and is quick in typical applications, 
where -/V gr id is the number of grid points for a given pa- 
rameter. Next, we map the corners of a hypercube to 
these neighboring points, 2 N in total for an AD hyper- 
cube. 

The weight of the theoretical value of a grid parameter 
point to that in the targeted parameter point is illus- 
trated in Figure El in the 2D case. The weight of point 
(rro, yo) to the targeted point is the area of the stippled 
square, \5xq ■ 6yo\. When the point (x,y) coincides with 
(#2, 2/2), Sxq -5yo = 0, or the weight on the targeted point 
value due to point (xq, yo) is zero. The final interpolated 
value at the targeted point is calculated by summing up 
contributions from all of the points 

f(x,y)=f(x ,y ) ■ \Sx Q ■ 6y Q \ 
+f(xi,yi) ■ \Sxi ■ Sy x \ 
+f(x2,V2) ■ \Sx 2 ■ Sy 2 \ 
+f(x 3 ,y 3 )-\Sx 3 -Sy 3 \ . (2) 

For a square of unit length for each side, we can re- write 
\5x ■ Sy \ to be (1 - \x- x \) ■ (1 - \y-yo\), similarly for 
other corners. 

To extend to the AD parameter space, the interpo- 
lation formula for the value of interest at the targeted 
point x, f(x), is 



2 N 

E 



j(z l )w\x) 



(3) 



where i denotes the sum over contributions from all of 
the neighboring points, z % , enclosing the targeted point, 
x, and the weight of each neighboring point is 



w\x) 



N 
3=1 



(4) 



in which j denotes a given parameter axis in the iVD 
parameter space. For a hypercube of length unity 
for each side, w l is equal to unity. This for- 
mula is an exten s ion of linear interpolation in the ID 
case by iBurnettl |l986, his Eqn. 5.2), and is shown 
to be applicable to computations in fluid dynamics 
(M urman. Aftosmis and Nemedl2004l ). 

The Eqns. [3] and 3] assume normalized coordinates, as 
such the range of each coordinate value lies within [0, 1]. 
Since the neighboring points z 1 are represented by the 
corners of a hypercube, no extra normalization step is 
required for those. The target point is normalized ac- 
cording to 

x j = T ~7 j (5) 

Pj -Pj 

where pj is the actual parameter value (e.g., age = 3 Gyr) 
at point Xj , and pj and p~j are that of the pair of neigh- 
boring points enclosing the targeted point, for the j-th 
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parameter axis (e.g., age = 1 Gyr and 4.2 Gyr, respec- 
tively). To derive the flux densities of a targeted spec- 
trum, Eqns. [3] and [4] are applied on the wavelength-to- 
wavelength basis. If the uncertainty in the flux density of 
the theoretical spectra is also available, Eqn. O together 
with usual error propagation formulae, can be used to 
obtained the uncertainty in the interpolated spectrum. 

The accuracy of the interpolated value at an intermedi- 
ate parameter point in-between the initial grid points de- 
pends on how good the linearity assumption is. This sim- 
plifying assumption is made only locally, i.e., in-between 
two grid points in the initial parameter grid, for a given 
parameter axis. There is no requirement that the linear- 
ity has to stand in the global range of a parameter. The 
reason is that any non-linear change of the theoretical 
value with parameter globally can be taken into account, 
e.g., by adopting a higher sampling fraction of the initial 
grid points in the related parameter regions. If indeed 
the spectral features depend highly non-linearly with pa- 
rameters that happen to occur at unknown parameter 
amplitude(s), a preliminary step can be taken to locate 
the concerned parameter amplit udes. The recently intro- 
duced locally linear embedding (jVanderplas fc Connolly] 
2009) appear to be a promising approach for this pur- 
pose, for its ability to unfold a non-linear manifold (see 
their Figure 1). 

3. EFFECT OF GRID RESOLUTION ON PARAMETER 
ESTIMATES 

A fundamental aspect in parameter estimation is the 
choice of grid resolution. A finer parameter grid is ex- 
pected to give higher accuracy in the estimates, whereas 
a smaller one requires less computational time. We use 
here the on-the-ny ability of our model generation ap- 
proach to investigate the effect of the targeted grid res- 
olution on the parameter estimates. This study requires 
multiple sets of models of increasingly large sample size. 

To construct an initial grid we use the 

iBruzual fc Chariot! (|2003l ) stellar population model 
(resolution 1 A per wavelength bin within the vacuum 
wavele ngth range 3450 - 8350 A), and the lCalzetti et all 
( 2000) intrinsi c dust model. T he initial stel lar mass func- 
tion i s that bv lChabrieil (|2003[ ) (similar to iKroupa et al.l 
119901 with a turnover below ~ 0.3 Mq) with lower 
and upper limits in stellar mass, 0.1 — 100 Mq. A 
wide range of stellar population and dust parameters 
is covered - the age of the oldest stars, age, takes the 
values: 1.0,4.2,7.3,10.5,13.7 Gyr, stellar metallicity, 
Z: 0.0004,0.004,0.008,0.02,0.05, the e-folding time 
of an exponentially decreasing star formation history, 
r e : 1.0,4.5,8.0,11.5,15.0 Gyr and the color excess, 
E(B - V): 0.0,0.2,0.4,0.6,0.8 mag. The resultant 
(4D) initial model grid is consist of 625 dust-attenuated 
integrated stellar spectra. The initial grid is fixed 
throughout the analysis. 

An example interpolated spectrum is shown in Fig- 
ure [3] Its age, Z, r e and E(B — V) are respectively 
8 Gyr, 0.02, 9 Gyr, 0.3 mag. The difference spectrum, 
interpolated - true theoretical spectra, is also shown. 
The true theoretical spe ctrum is referred to that being 
output directly from the IBruzual fc Chariot] (I2003D pro- 
gram and post dust-attenuated with the iCalzetti et al.l 
(2000) model. The flux density error is driven by both 
the parameter resolution in the initial grid and the 



validity of the linearity assumption within a single 
bin of the initial grid. In the interpolated spectrum, 
the error is about 2%, assigned equally to each and 
every wavelength bin to give a reduced % 2 between 
the two spectra close to unity (equals 1.07 in this 
case). Extending to the full parameter space, we gen- 
erate 320 spectra at age — 3.0,6.0,8.0,12.0 Gyr, 
Z = 0.0004, 0.004, 0.008, 0.02, 0.05 3 , r e 
3.0,6.0,9.0,14.0 Gyr and E(B — V) = 

0.1,0.3,0.5,0.7 mag using the above-mentioned ini- 
tial grid, which are compared with the true theoretical 
spectra. Figure |4] shows the distribution of the flux 
density error in the interpolated spectra, showing an 
average of 3.5%, a minimum of 0.4%, and a maximum 
of 11.3%; with a resultant average \ 2 °f 1-0- The 
amplitude of the flux density error is comparable with, 
e.g., the spectrophotometric uncertainty from the Sloan 
Digital Sky Survey Data Rele ase 6, 7% in the observed 
frame wavelength of 3800 A ( Adclman-McC arthv et aLl 
2008). The threshold in the flux density error, at which 
the interpolated spectrum is still useful, would depend 
on individual applications. We plan to extend this type 
of comparison t o other models . One model to consider, 
e.g., is that by iKotulla et al.l (|2009h . who provide an 
interactive web interface for downloading theoretical 
spectra of integrated stellar systems which are defined 
by multiple parameters. 

We next construct a mock galaxy sample for the pa- 
rameter estimation, in which the parameters underlying 
each spectrum are known. Based on the initial grid, a 
random sample of 100 theoretical spectra are generated 
through our approach. Random Gaussian noise is added 
to the flux density in each wavelength bin of each spec- 
trum, fixed at signal-to-noise (S/N) of 20. Bayesian 
parameter estimation is then performed on the mock 
sampl e, using the approach described by iKavirai et al.l 
(2007), where all of the 4 model parameters are esti- 
mated simultaneously. The parameter estimations are 
carried out using several different samples of theoretical 
spectra, defined from high to low parameter resolution, 
as shown in Table [1] These spectra are also generated 
using the model generation approach described in this 
work, and are based on the same initial grid. As such, we 
are probing the effect on the parameter estimates due to 
the increasing targeted grid resolution only, rather than 
that due to model imperfection in terms of describing 
the real observed spectra, and so o n. Except for the stel - 
lar metallicity, for which we follow iGallazzi et all {2005) 
and adopt logarithmic bins (i.e., evenly sized in log 10 Z), 
linear bins (i.e., evenly sized bins) are used for all of the 
other parameters (age, T e and E(B — V)) in the targeted 
grid. 

The dependence of the root mean square error 
(RMSE), the square root of the variance and the bias 
(the average difference between the estimate and the true 
parameter value) of the estimated parameter are shown 
in Figure [5] for each of the parameters. Mathematically, 
the RMSE is the square root of the mean square error 

3 The stellar metallicity values chosen here are not different from 
those in our initial grid , because other metallicity values are not 
output directly from the Bruzual & Chariot (2003) program. So in 
this example we are in fact considering a 3D interpolation in a 4D 
hypercube. 
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Fig. 3. — An example i nterpolated spectrum (vacu um wavelength) by using a 4D hypercube, shown in red. The true theoretical spectrum 
from the dust-attenuated Bruzual & Chariot (2003) model is plotted for comparison, in blue. The difference spectrum (interpolated - true 
theoretical) is shown in green. Some Balmer absorption lines are marked. The solar luminosity Lq is 3.826 X 10 33 ergs s — 1 . The flux 
density error in the interpolated spectrum is about 2%, assigned to give a reduced \ 2 close to unity. 
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Tested parameter grids. 
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Fig. 4. — The distribution of the error in the flux density of the 
interpolated spectra over a 4D grid, wit h respect to the tru e theo- 
retical spectra from the dust-attenuated Bruzual & Chariot (2003) 
model. The vertical dashed line indicates the average error, 3.5%. 

(MSE). The MSE is the sum of the variance and the 
bias squared, so it can be interpreted as t he sum of th e 
squares of statistical and systematic errors (|CowanlH9 98). 
The MSE is found to decrease as the resolution of the 
targeted parameter grid increases, an indication of an 
improved determination of the posterior probability for 
each parameter of the mock galaxies, in this case, due 
to the increase in the number of parameter grid points 
being sampled. We also see that the estimates are fairly 
unbiased when the number of theoretical spectra > 10 4 
for our parameter ranges, where the corresponding grid 



Targeted grid size 


Reduction factor 


Number of spectra 


14 x 36 x 15 x 17 


1 


128520 


10 x 24 x 10 x 12 


1.5 


28800 


7 x 18 x 8 x 9 


2 


9072 


4x9x4x5 


4 


720 



There are 4 parameters (age, Z, r e , E(B — V)) constituting each 
of the grids, in that order the grid size is given. The range of a 
parameter is the same for all of the grids, given in [[3] The reduc- 
tion factor is the approximate fractional decrease of grid resolution 
relative to the case of highest resolution, and is set to be uniform 
in each parameter. Except in Z for which logarithmic bins are 
adopted, linear bins are used in all of the other parameters. 

resolution is listed in Table [1] On the other hand, the 
variance of the estimates at our large-sam ple limit should 
be lim ited by the data S/N. For example. IGallazzi et alj 
(2005) compared (their Figure 3) both the shape and 
the width of the posterior probability in stellar metal- 
licity between the high- and law-S/N spe ctra from the 
Sloan Digital Sky Survey (|York et al.ll2000f ) , in that lower 
S/N spectrum gives wider distrib ution in the posterior 
proba bility of parameters. Earlier [Gil de Paz fc Madord 
(2002) also came to similar conclusion, through parame- 
ter recovery test upon synthetic photometry of integrated 
stellar populations. 

Another interesting aspect of Figure [5] is the simul- 
taneous improvement in all of the parameter estimates 
upon the increase in the targeted parameter grid resolu- 
tion. This result manifests the presence of degeneracies 
among all of the involved parameters - age, Z , r e and 
E(B — V). It is well-known that age-metallicity degen- 
eracy may hinder the studies of astrophysical objects by 
using only their broadband colors. Even integrated spec- 
tra are not immune from this degeneracy - an integrated 
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Fig. 5. — The dependence of the RMSE (filled black circles), the 
square root of the variance (filled gray circles) and the bias (empty 
circles) of the parameter estimates on the number of theoretical 
spectra, or the resolution of the targeted parameter grid (see Ta- 
ble [TJ. The solid horizontal lines mark the zero value. The dotted 
lines are for eye-guiding only. 

stellar spectrum with intermediate to old age (1.5 Gyr 
and up) may look the same by tripling the age or reduc- 
ing t he stellar metallicity by a factor of two (jWorthev! 
1994). To determine age, one therefore needs to know 
Z , and vice versa, meaning both parameters have to be 
determined simultaneously. Figure [5] demonstrates that 
degeneracies exist among all of the considered parame- 
ters, in the way that the improvements on these estimates 
go hand-in-hand. 

4. SUMMARY 

In applications where a large (sample size > 10 5 ), non- 
analytical model is needed, the hypercube representation 



of the parameter space, together with multi-linear inter- 
polation for deriving the theoretical value at an interme- 
diate grid point, are shown here to be feasible to improve 
the flexibility in handling the models. Using stellar pop- 
ulation synthesis of galaxies as a pilot study, we found 
that the increase in targeted parameter grid resolution 
improves the parameter estimates in terms of the mean 
square error. This example is given by incorporating the 
Bayesian parameter estimation. Naturally, our model 
generating approach can be combined with other param- 
eter estimation techniques, such as the Markov Chain 
Monte Carlo. 

Models of large sample size can be applied to studies in 
various areas and disciplines. In the context of spectral 
analyses, a large grid of theoretical spectra is expected 
to be applicable to mock catalog construction, study- 
ing the relationship between spectral features and model 
parameters, spectral fitting, and stellar population syn- 
thesis of galaxies. Our approach is particularly suitable 
for studies that are based on multi-dimensional param- 
eter space, and are focused on investigating differences 
among results obtained through different models. 
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